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ABSTRACT 

We have derived the first, fuUy-cosmological, similarity solutions for cold dark mat- 
ter (CDM) halo formation in the presence of nongravitational collisionality (i.e. elastic 
scattering) , which provides an analytical theory of the effect of the self-interacting dark 
matter (SIDM) hypothesis on halo density profiles. Collisions transport heat inward 
which flattens the central cusp of the CDM density profile to produce a constant- 
density core, while continuous infall pumps energy into the halo to stabilize the core 
against gravothermal catastrophe. This is contrary to previous analyses based upon 
isolated haloes, which predict core collapse within a Hubble time. These solutions im- 
prove upon earlier attempts to model the formation and evolution of SIDM haloes, offer 
deeper insight than existing N-body experiments, and yield a more precise determina- 
tion of the dependence of halo density profile on the value of the CDM self-interaction 
cross section. Different solutions arise for different values of the dimensionless colli- 
sionality parameter Q = apbrvir oc r^ir/Xmip, where a is the scattering cross section 
per unit mass, pb is the cosmic mean matter density, rvir is halo virial radius and Amfp 
is the collision mean free path. The maximum flattening of central density occurs for 
an intermediate value of Q, Qth, at which the halo is maximally relaxed to isother- 
mality. The density profiles with constant-density cores preferred by dwarf and low 
surface brightness galaxy (LSB) rotation curves are best fit by the maximally-flattened 
{Q = Qth) solution. If we assume that dwarfs and LSB galaxies formed at their typical 
collapse epoch in ACDM, then the value of a which makes Q = Qth is cr ~ 200 cm^ g~^ , 
much higher than previous estimates, a ~ [0.5 — 5] cm^ g""'^, based on N-body experi- 
ments. If a is independent of collision velocity, then the same value a ~ 200cm^g~^ 
would make Q > Qth for clusters, which typically formed only recently, resulting in 
relatively less flattening of their central density profile and a smaller core. 
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1 INTRODUCTION 

Tlie cold dark matter (CDM) model provides a successful 
framework for understanding the formation and evolution 
of structure in our universe. According to this model, grav- 
ity amplifies primordial density fluctuations, and structure 
forms by hierarchical clustering: small objects form first, 
later merging to form larger objects. The most promising 
candidate for CDM is the weakly interacting massive par- 
ticle (WIMP). In this picture, the microscopic interaction 
between CDM particles is negligible (coUisionless) , and they 
interact only by gravity. However, this assumption has not 
been fully verified and it is fair to say that the microscopic 
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nature of CDM is still unknown. It is important, therefore, to 
explore the consequences of varying this underlying assump- 
tion about CDM in the hope that astronomical observations 
can be used to place meaningful constraints. 

The possible variation of the microscopic nature of cold 
dark matter is closely linked to the problems of the CDM 
model. Despite its success, the CDM model has several prob- 
lems w hich exist mo stly in the small scale regime (see, for ex- 
ample, [Moot^IQ^)- Among these problems, much attention 
has been focused on the N-body simulation results for the in- 
ner density slope, since the observed rotation curves of dark- 
matter dominated dwarf and low surface brightness (LSB) 
disk galaxies tend to favour mass profiles with a flat-density 
core unlik e the singular profiles of the CDM N-body simula - 
tions (e.g. iFlores fc Primac^ll994l: iMarchesini et al.ll2002ll . 
The latter are generally characterized by an empirical fitting 



© 0000 RAS 



2 K. Ahn & P. R. Shapiro 



formula for the spherically averaged density profiles in those 
N-b ody results, either the Navarro- F renk- White (NFW) pro- 
file (|Navarro. Frenk fc Whitel |l gg7|). for which o oc as 
r 0. or the Moore profile iMoore et al.lll999^ . for which 
p oc r^^'^ instead^. It was controversial whether the ob- 
served d ata was resolved well enough to indicate a soft core 
(e.g. see Ivan den Bosch fc Swatersll200ll for possible beam 
smearing effect on the data), but observations have built up 
which favour the soft core even after elimin ating the beam 
smearing effect (e.g. iMarchesini et afll2002l and references 
therein). After much more work on the N-body simulations 
of CDM, the discrepancy between these data and the nu- 
merical halo profiles remains si gnificant (see, for instance, 
iDiemand. Moore fc Stadelll2004l) 

This apparent discrepancy between the N-body sim- 
ulation results and observations of dwarf galaxy rotation 
curves has raised a question about the nature of dark 
matter, including the assumption that CDM is collision- 
less. People have suggested solutions to this discrepancy 
which either preserves the coUisionless nature of the dark 
matter or else adopts a new picture. In the former cate- 
gory are explanations which attribute the discrepancy to 
astr ophysical processes beyond the pure N-body dynamics 
fe.fi.lEl-Zant. Shlosman fc Hoffmanl200ll : fWeinberg fc Katj 
I2OO3) or to a primordial power spe ctrum tilted away 
from the Harrison-Zel'dovich shape fe.g. lZentner fc Bullockl 
|20o3). In the latter category, the proposal of self- interacting 
dark matter (SIDM) bv iSoergel fc Stcinhardt (200(1) has re- 
ceived a lot of attention. In this picture, microscopic inter- 
action between dark matter particles is non-negligible and 
can affect the dynamics of halo formation. Since the actual 
identity and microscopic nature of CDM is still unknown, it 
is important to explore the consequences of such hypotheses 
in the hope that astronomical observations can be used to 
place meaningful constraints. 

The trunca ted isothermal sphere ( TIS) halo model 
developed by IShapiro. Iliev fc Raeal lll999l . also see 
llliev fc Shapirdl200lll provides a simple physical clue about 
the existence of such soft cores in haloes of cosmological ori- 
gin which otherwise closely resemble the CDM haloes found 
by N-body simulation in all respects except the presence of 
the inner cusp. The model is based upon the assumption 
that haloes form from the collapse and virialization of "top- 
hat" density perturbations and are spherical, isotropic, 
and isothermal. This leads to a unique, nonsingular TIS, 
a minimum-energy solution of the Lane-Emden equation. 
The resulting density profile is nearly indistinguishable 
from that deduced from the observed ro tation cu r ves o f 
dwarf and LSB disk galaxies, as fit bv [ Burkerd (Il995l) 
with a profile involving a soft core jlliev fc Shapiro! I2OOII) . 
The assumption of isothermality which underlies the TIS 

1 Recently, iHavashi et"all <2004D and iNavarro et"ai] i2004D re- 
ported that the logarithmic slope of the spherically averaged den- 
sity profile of dark matter haloes in their ACDM simulations, 
— dlnp/dlnr, decreases monotonically towards the centre. This 
change of slope, however, does not seem to fully account fo r the 
observed rotation curves. iDiemand. Moore fc Stadell ^2004^ con- 
clude, by surveying several simulation results by different groups, 
that neither the typical NFW profile nor the Moore profile pro- 
vides a good fit, but one ne e ds a little more complex form such 
as the one bv lNavarro et all i2004h . 



halo model is the primary reason that the central cusp of 
the CDM N-body haloes is replaced by a soft core. The 
N-body haloes are also approximately isothermal, but there 
is typically a small temperature dip near the centre. This 
suggests that if some mechanism existed to transport heat 
inward so as to make CDM haloes more isothermal, they 
might exhibit soft cores while maintaining the same basic 
halo structure at larger radii already found for CDM haloes. 
The coUisionality of SIDM serves as such a mechanism: 
elastic collisions transport heat inward, which flattens 
the central cusp of the CDM density profile to produce a 
constant-density core. 

The problem of SIDM halo formation has so far been 
studied primarily by numerical N-body experiments. Some 
of the first attempts to calculate the effect of the elastic 
scattering of SIDM on halo structure formation involved 
isolated haloes which were assumed initially to follow the 
equilibrium profiles (e.g. NFW profile) found in collision- 
less N-body simulations of standard CDM ( Burkcrt 20o3; 
Kochanck fc White 2001; both start from the Hernquist 
profile ( Hernquist , 1990 | ) which is a go o d app roximation to 
the NFW profile). iKochanek fc Whitd (l200ll) . for instance, 
found that SIDM haloes can form flat density cores within a 
relaxation time as expected, but also that the lifetime of such 
flat cores is only a few relaxation times. They concluded, 
therefore, that most galactic haloes would have undergone 
core collapse. 

The effect of SID M coUisionality on halo structure 
has also been studied bv lBalberg. Shapiro fc Inaeakil ||2002| . 
BSI hereafter) by solving ID, quasi-static fluid equations. 
These authors also considered isolated haloes like those in 
the N-body experiments mentioned above, adopting ID, 
spherical symmetry, with non-cosmological boundary con- 
ditions. They treated the dynamics of SIDM by a fluid ap- 
proximation developed previously in the study of stellar dy- 
namics, derived from the Boltzmann equation, in which they 
modifled the heat conduction term to handle the elastic scat- 
tering of the SIDM particle-particle interaction. They solved 
the spherically-symmetric, virialized "gravothermal fluid" 
equations of |Lvnden-Bell fc^Eggletojl (1980), which include 
mass conservation, hydrostatic equilibrium, an equation for 
heat conduction, and the flrst law of thermodynamics. Ac- 
cording to these equations, the halo is time-dependent be- 
cause heat conduction causes it to evolve through a contin- 
uous sequence of hydrostatic equilibria. An analytical self- 
similar solution to these equations was found in the limit of 
large scattering mean free path - the limit where the mean 
free path is much larger than the size of the halo - following 
the derivation of iLvnden-Bell fc Eggletoiil il98Ci) for globu- 
lar clusters. The BSI solution shows that secular evolution 
always takes a conflguration in the long mean free path limit 
and drives it into the short mean free path regime. In order 
to track the evolution into the short mean free path regime, 
as well, BSI then used their similarity solution as the initial 
condition for a numerical solution of the same gravother- 
mal fluid equations, but for finite scattering cross-section. 
This approach made it possible to follow core collapse to 
a much more advanced stage than the N-body experiments 
could. From this, they concluded that the ultimate core col- 
lapse time was much larger than the relaxation time in the 
core, long enough even to exceed a Hubble time in some 
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cases^. Their estimated core collapse t ime is ^c o u — 2 90 
which contradicts the resu lt found by iBurkerd i200Cfl and 
iKochanek fc Whitel l|200ir i that the core collapse time was 
only a few relaxation times. 

This apparent discrepancy in core collapse timescale 
between t he study by BSI and the numerical N-b ody ex- 
periments ljBurkertll2nnnl: iKochanek fc WhitJ20nl^ may be 
attributed to the fact that their adopted initial conditions 
were different. In BSI, as mentioned, the initial condition 
was tuned to be in the extremely long mean free path limit, 
^mip/H 3> 1, where A^fp is the mean free path and H 
is the gravitational scale height or roughly the halo size. 
This occurs when the system is either dilute enough or the 
scattering cross section (per unit mass) g is small enough, 
since A^fp oc l/{ap). According to BSI, most of the col- 
lapse time is spent to reach the condition Amfp/i? — 1, 
and the halo density profile always has a fiat core. By 
contrast, Kochanck & White (2001) started with a cuspy 
profile with parameters which corresponded to a condition 
^mtp/H ~ 0.1—3.0. In this case, SIDM halo cores can quickly 
reach \^{p/H ~ 1 or they are already in the short mean free 
path limit, which then requires only a few relaxation times 
for core collapse. 

In what follows, we will cover the whole range of a 
and show that the observed dwarf-galaxy rotation curves are 
best-fit when the SIDM interaction has its maximal effect, 
which occurs when Xraf-o /H ~ 1, so the r egime of greatest 
interest may be that of IKochanek fc Whit e (2001J. If so, 
then most isolated haloes would, indeed, suffer core collapse 
within a Hubble ti me. However, the shared limitation of 
the analyses of BSI. lBurkerl J2000h and IKochanek fc Whitj 
l|200ll^ . that of non-cosmological boundary conditions, is a 
severe one. Cosmological infall may inhibit core collapse. If 
cosmological infall can delay core collapse substantially, pre- 
vious estimates based on isolated halo models would change 
by shifting the time of the onset of core collapse until cos- 
mological infall becomes negligible. 

The effect of cosmological boundary conditions on 
SIDM halo formation has been studied numerically by cos- 
mological N-body simulations, in which Gaussian random 
noise initial conditions for CDM were incorporated. Early 
work along these lines attempted to derive the maximal ef- 
fect of coUisionality by adopting the fully c oUisional limit 
which corresponds to ordi nary gas dynamics jYoshida et alJ 
l2000al: iMoore et alll200tt . The surprising resuh they re- 
ported was that simulations yielded density profiles with 
central cusps even steeper - with logarithmic slope close 
to —2 - than those in coUisionless N-body simulations. 
Subsequent cosmological N-body experiments which treated 
the SIDM elastic scattering in more detail, however, re- 
ported that halo density profiles were fiatter than those 
of either fully coUisional or purely coUi s ionless simulation s 
jYoshida et al.ll2000bl: iDave et alJl200ll: IColin et al.ll2002l) . 
They found that values of a in the range a ~ [0.1 — 5] cm'^g"^ 
(the range of preferred a varies among these works, but 
within less than an order of magnitude) produced SIDM 
haloes with sufficient profile fiattening to account for dwarf 



In a follow-up paper, the authors showed how this process 
could lead to the formation of seeds for super-massive black holes 
jBalberg fc Shapiro 2002). 



galaxy rotation curves, which did not suffer from the rapid 
core collapse identi fied in the earlier N-body experiments 
for isolated haloes. iDave et alJ l)200 J) speculated that this 
might be the result of cosmological infall, which was ab- 
sent from the calcula tions of isolated haloes. In addition, 
lYoshida et al.l l^2000b^ found that SIDM collisions in cluster- 
sized haloes would be more frequent than in dwarf-sized 
haloes, thus producing relatively larger cores in clusters, 
while observations tend to show that dwarfs and LSBs show 
relatively larger cores than clusters. This led them to sug- 
gest that the scattering cross section be velocity-dependent 
(a oc 1/v) so that more massive haloes would have a smaller 
degree of density profile flattening. Later, jColin et al. (200 J 
tested this hypothesis and confirmed that it could match ob- 
served core sizes from dwarfs to galaxy clusters. 

Further study of the formation and evolution of SIDM 
haloes is warranted to resolve the issues raised by previous 
work and put the subject on a firmer theoretical footing. 
On the one hand, the numerical N-body simulations and 
ID semi-analytical treatment mentioned above of isolated 
haloes with non-cosmological boundary conditions are un- 
able to address the important effects of cosmological infall. 
On the other hand, the fully cosmological N-body simula- 
tions which have been performed of haloes that arise dur- 
ing large-scale structure formation in the SIDM model have 
so far been limited by numerical resolution and dynamic 
range. As a result, simulation results published to date do 
not attempt a detailed enough comparison with observed 
galactic rotation curves to determine if SIDM haloes can re- 
ally match them or to give a very reliable constraint on the 
SIDM cross section. Such simulations do not afford enough 
insight into the underlying dynamical processes which gov- 
ern the halo structure, either. Finally, the wide range of 
typical collapse epochs expected in a CDM universe for the 
wide range of halo masses extending from dwarf galaxies to 
clusters suggests that the effect of the SIDM interaction may 
be halo-mass-dependent; this dependence has not yet been 
adequately explored. 

Towards this end, we have derived a fully cosmological 
model for the origin and evolution of CDM haloes in the 
presence of nongravitational coUisionality (i.e. elastic scat- 
tering). We have combined the fluid approximation of BSI 
with the spherical infall model for cosmological perturbation 
growth to yield fully time-dependent, detailed similarity so- 
lutions for SIDM haloes, for arbitrary degree of coUisionality. 
We shall apply these solutions to test the hypothesis that 
cosmological infall retards the core collapse of SIDM haloes, 
and compare the predicted SIDM halo proflles with the mass 
profiles inferred from dwarf galaxy rotation curves. This will 
enable us to place much better quantitative constraints on 
the SIDM cross section and better assess the validity of the 
SIDM model. 

S ubsequent to the original suggestion by 

ISpergel fc Steinhardd 120001 and the exploration de- 
scribed above of the halo structure which results, related 
work has focused on constraining the SIDM hypothesis 
by its implications for other astrophysical phenomena or 
attributing density flattening to more complicated CDM 
dynamics or gas dynamical feedback effects within the stan- 
dard CDM picture. There seem to be s trong observational 
constraints on the possible range of a. ICnedin fc Ostrikerl 
(2001) ruled out a range of a, cr = [0.3 — 10'']cm^g~^, 
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based on their calculation of the evaporation time of the 
dark matter h aloes of elliptical galaxies in the clusters. 
iNataraian et al. (2002) rule out all values of ct > 42cm^g~^ 
by comparing the predicted truncation radii of SIDM 
haloes inside clusters by ram-pressure stripping to those 
of observed haloes which they obtain using cluster gravi- 
tational lensing observations. iHennawi fc Ostrikeil (l2002l) 
conclude that if cr 2> 0.02 cm^g~^, the supermassive black 
holes in the centres of galactic haloes would be more 
massive than observed. Along the line of explaining 
flattened cores within the standa rd CDM picture (se e 
iPrim ackl l20 0j for recent review). lEl-Zant et alj J2004l) . 
as in El-Z ant. Shlosman fc HoffmanI (|20mF ! investigate 
the effect of dynamical friction by clumpy substructure 
associated with baryonic dissipation and show th at a 
flat core can be generated. IWeinberg fc Kat j i2002h and 
iHoUev-Bockelmann. Weinberg ^^<£^all20o3r claim that 
soft cores can be induced by the presence of a CDM 
bar structure. iHavashi et al.l ll2004) argue that the gas 
rotation speed in a galactic disk may be different from 
the dark matter circular velocity, which would, therefore, 
be wrongfully interpreted as implying the existence of 
soft cores. Dynamical feedback ef fects from supernova 
explo sions is another possibility iNavarro. Eke fc FrenkI 
Il997l) . Observationally, related to the missing satellite 
problem, the gravitational lensing flux anomaly seems to 
require clumpy, da rk substructures fMetcalf & Zhao" 2002"; 
Dalai fc Kochaneh 2002.; .Keeton. Gaudi fc Fetters 2003,; 



Mao et al.ll2004l) . even though no firm conclusion has been 



established yet. 

Nevertheless, we believe that SIDM is still a viable can- 
didate for dark matter. As we shall discuss below in 
the previous analyses restricting cr are subject to significant 
caveats. It is meaningful to study this subject because it can 
shed light on the nature of dark matter whose origin we do 
not know yet, and simply because CDM problems are far 
from reaching a firm conclusion. Also, the research effort on 
SIDM has not yet reached the level of that on coUisionless 
CDM. 

We have previously reporte d the resu lts of this paper 
in a summarized form elsewhere. |Ahn fc Shapiro 120Q3a b) 
briefly described the result from H5.5.2I Sh apiro et alJ (I20M1 
summarized the results in light of the current status of CDM 
research. However, we warn the reader that these works 
should be considered as guidelines to the detailed, complete 
picture of this paper. 

The plan of this paper is as follows. In f|5| we review 
the previous work on self-similar gravitational collapse in 
the spherical infall model. In we show how the fluid ap- 
proximation derived from the Boltzmann equation can be 
generalized to include a conduction term which accounts 
for nongravitational coUisionality in an otherwise collision- 
less dark matter gas. The assumptions which underlie this 
fluid approximation will be justified by comparison with the 
results of cosmological N-body simulations of CDM haloes 
and with detailed similarity solutions for halo formation by 
gravitational collapse of coUisionless dark matter. In we 
show that gravitational collapse by spherical infall can be 
self-similar even in the presence of coUisionality, and that 
the specific condition required for self-similarity is relevant 
to galactic haloes. We then describe the basic equations for 
the dynamics of SIDM haloes. In we present the simi- 



larity solution and its application to test the SIDM hypoth- 
esis, including a comparison with previous N-body results 
and observed galactic rotation curves. Our conclusions and 
discussion are in 



2 SELF-SIMILAR GRAVITATIONAL 
COLLAPSE 

2.1 Previous analytical models for halo formation: 
the spherical infall model 

Analytical approximations have been developed to model 
the formation of haloes by the ID growth of spherical cos- 
mological density perturbations, in either a coUisionless gas 
or a fluid. We shall need to refer to some of these solutions 
to justify our fluid approximation in S|3 We shall also build 
upon these earlier solutions in deriving new ones which add 
the effects of coUisionality in >13.2l andl^ It is necessary for us 
to b egin, therefo re, by briefly recounting this earlier work. 

iGunn fc Go tt (197 J) first presented the concept of the 
so called "secondary infall model (SIM)". This SIM refers to 
the effect of the addition of a point mass to a uniform, ex- 
panding Friedmann-Robertson- Walker universe as a pertur- 
bation which causes the surrounding spherical shells to de- 
celerate relative to the background universe, until they reach 
a radius of maximum expansion and re-collapse. Subse- 
quent work generalized this approach to include spherically- 
symmetric initial perturbations for which the overdensity 
profile depends upon radius or mass as a scale-free power- 
law. Along these lines, nFillmorc fc Goldroich (1984,, FG here- 
after) studied the dynamics of coUisionless CDM haloes 
using a self-similar model, adopting a scale-free initial 
overdensity param etrized by its shape; e in equation 10. 
'Bcrtschin geJ ^1985^ studied a special case of FG but also ex- 
tended the analysis to a coUisional fiuid. . Hoffman fc ShahamI 
(1985, HS hereafter) showed that a power-law power spec- 
trum would indeed generate a scale- free initial condition, 
such as was adopted by FG. They then argued that the re- 
sulting nonlinear structure would be described by a power- 
law profile determined by the shape of the power spectrum. 

Previously mentioned works adopted a rather un- 
realistic condition for the coUisionless case, that of 
purely radial motion. N-body simulations of CDM halo 
formation find that the virialized region tends toward 
isotropic random velocities. Some attempts to incorporate 
tangential velocities within the frame work of spheri- 
cal symmetry hav e also be en made by Rvden fc GunnI 
(1984)^ Avila-R ces. Firman i fc Hernandez (1998) and 
iHozumi. Burkert fc Fuiiw ara (2000). Alo ng these lines, one 
may a ls o refer to work by .Kull L19 99I). iLokas fc HoffmanI 
J2000l) . iDel Popolo et al.! J2000l) . iHiotelisI J2002|) and 
references therein. 

The fluid approximation of coUisionless CDM halo 
forma tion has emerged recently. iTevssier. Chieze fc Aliml 
(Il997f) showed that one could use fluid-like conser- 
vation equations t o mimic the radial-only SIM s such 
as the FG model. ISubramanian. Cen fc Ostrikerl (l200(Jl 
extended this analysis to include tangential motion. 
In these pictures, one could expect, for instan c e, ac- 
cret ion shock structure. iTexssiglj Chieze fc AUmil (ll99iD 
and ISubramanian. Cen fc Ostrikenil200ol) therefore used a 
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"shock jump condition" which one would only expect from 
a purely collisional fluid. The fluid approximation has been 
used in the literature of stellar dynamics, but its applica- 
tion to CDM halo dynamics is rather new, and as will be 
described in it simplifles the description of dark matter 
dynamics substantially. 

Self-similar spherical infall models have also been used 
to stud y the effect of inc o rporating additional baryonic 
physic s. [Bertschinger ('1989*) , 'Owen. Weinberg fc VilhimserJ 
il998r . and Abad i. Bower fc Navarr o (2000) have studied 
the effect of gas cooling on galaxy formation using self- 
similar models. Because the system was tuned to maintain 
self-similarity in the presence of cooling, the cooling func- 
tion used is not perfectly physica l. Nevertheless, as shown by 
lAbadi. Bower fc Navarrd l)200nh . one can use these models 
at least to test hydrodynamic codes in the presence of cool- 
ing. Moreover, these models capture the generic behaviour 
of galactic dynamics in the presence of realistic cooling. 

The similarity model we shall develop in this paper is 
the first self-similar halo model which includes the effec- 
tive heat conduction resulting from SIDM coUisionality. In 
section JOandIO we describe how SIMs with scale-free 
linear perturbations arise naturally from the Gaussian ran- 
dom noise primordial density fluctuations, thus leading to 
the self-similar evolution of nonlinear structures. 



2.2 Halo formation from scale-free linear 
perturbations 

In the Einstein-de Sitter (EdS) background universe, an ini- 
tial linear perturbation whose mass profile is spherically 
symmetric and has a scale-free, power-law form 

5M 



oc M~ 



results in structure formation which is self-similar (EG). 
Each spherical mass shell around the centre expands until 
it reaches a maximum radius (turnaround radius rta), and 
re-collapses. For a given e, we have 

rta oc t^, 



where 

2 f3e + l- 
3 



(2) 



(3) 



Since there are no characteristic length or time scales for 
this problem other than the turn-around radius rta and the 
Hubble time t, the gravitational collapse which ensues from 
this scale-free initial condition must be self-similar as long as 
the background universe is Einstein-de Sitter, in the absence 
of physical processes which introduce additional scales (e.g. 
SIDM coUisionality) 

In general, if the unperturbed matter is a cold fluid, 
the infall which results from this perturbation is highly 
supersonic and is terminated by a strong accretion shock 
which thermalizes the kinetic energy of collapse. The ac- 
cretion shock radius is guaranteed by self-similarity to be 
a fixed fraction of rta(t) at all times. The mean density of 
the postshock region is, therefore, always a fixed multiple 
of the cosmic mean matter density. For most cases of inter- 
est, this postshock region is close to hydrostatic. For a col- 
lisionless g similar description applies as long as the in- 
falling matter initially had small (or zero) random motions. 
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Figure 1. Self-similar collisionless halo formation for e = 1: Com- 
parison of the skewless-fluid approximation to the exact colli- 
sionless Bertschinger solution. Solid lines represent the solution 
obtained from the fluid approximation in the radial direction, 
while dotted lines represent the collisionless Bertschinger solu- 
tion. Spikes in the density plot simply represent infinite values, 
corresponding to caustics, and therefore there is no physical sig- 
nificance in the height of these spikes. However, spikes in the 
velocity plot are finite and real. Note that solid lines do not rep- 
resent the 7 = 5/3 fluid Bertschinger solution. 



In that case, each mass shell collapses supersonically as a 
single stream until it encounters a region of shell-crossing 
and density caustics, which encompasses all previously col- 
lapsed (i.e. interior) mass shells. All collapsed mass shells 
inside this region oscillate about the centre. The radius of 
this region of shell-crossing, given by the outermost density 
caustic, is analogous to the shock radius in the fluid case. 

Results for the purely collisionless case w ere presented 
for arbitrary values of e by EG, and for e = 1 bv lBertschingeJ 
(, 1985) (where the latter included a fluid component, as well). 
Figures Q and 121 show the exact similarity solutions for the 
purely collisionless cases with e = 1 and e = 1/6, respec- 
tively. As we describe below in H2.3I these values roughly 
bracket the range relevant to cosmological haloes in a CDM 
universe. We will refer to these solutions again for compari- 
son in deriving our fluid approximation in ^ 



2.3 Halo formation from peaks of the Gaussian 
random noise primordial density fluctuations 

The theory of halo formation from peaks in the density 
field which result from Gaussian-random-noise initial den- 
sity fluctuations draws an interesting connection between 
the average density profile around these peaks and the shape 
of the fluctuation power spectrum. According to HS, local 
maxima of Gaussian random fluctuations in the density can 
serve as the progenitors of cosmological structures. They 
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10-3 lQ-2 10-' 1 10-3 10-2 10-1 1 



A A 

Figure 2. Same as Figure □] but £ = 1/6. Note again that the 
solid line was not generated from the 7 = 5/3 fluid approximation, 
but rather from the radial-only fluid approximation described by 
equations 1^ - E^ . 




4 6 8 10 12 14 



log,„ M (h-i M,,) 

Figure 3. Effective index of the power spectrum (P{k) oc 
vs. halo mass for the ACDM universe. Appendix describes how 
ricff is calculated. The solid line is derived from the HS approach, 
which we use in this paper. The dotted line is derived from cruj. 

for M in the range from W^Mq to 10"Mq). For haloes in 
the cluster mass range, M ~ IO^^Mq, n ~ —1.5. 



show that rare density peaks {v > 3, where f corresponds 
to uaM peak) have a simple power-law profile^ 

Ao(r)ocr-("+^), (4) 

where Ao(r) is the accumulated overdensity inside radius 
r, and n is the effective index of the power spectrum P{k) 
approximated as a power-law P{k) oc k" at wavenumber 
k which corresponds to the halo mass as described in Ap- 
pendix**. The overdensity Ao(r) is equivalent to the frac- 
tional mass perturbation 5M/M inside radius r, 

Ao(r-) = SM/M (xM~^. (5) 

From equations (0 and @, we deduce that the power-law 
power spectrum naturally generates a scale-free initial con- 
dition with 

e = in + 3)/3. (6) 

According to this model, as described in Appendix, 
haloes of a given mass M originate from density pertur- 
bations given by equation @ with n determined by the pri- 
mordial power spectrum after it is transferred according to 
the parameters of the background universe and the nature of 
the dark matter. We plot this effective n as a function of halo 
mass in Figure|21for the current ACDM universe. The value 
of n ~ —2.5 is a reasonable approximation for galactic haloes 
(i.e. n ~ -2.5±0.1 for M ~ IO^^^Mq, while n ~ -2.5±0.2 

^ iBardeen et al.l ll986l) also get a similar result: local density 
maxima have a triaxial profile, but as f increases it becomes more 
and more spherical with a profile converging to equation Bl . 
* The average linear overdensity profile, equation jl}, holds for 
any value of v. For small u, however, random dispersion around 
this average profile becomes substantial, limiting the generality 
of equation IH . 



3 FLUID APPROXIMATION OF DARK 
MATTER DYNAMICS 

We show that fluid conservation equations for a gas with 
adiabatic index 7 = 5/3 are a good approximation to the dy- 
namics of both CDM and SIDM haloes. This appro ach has 
been used in the literature of stellar dynamics (e.g. iLarsonl 
Il970l: iLvnden-BeU fc EggletonI Il98d : iBettwieseJ for 
the study of the gravothermal catastrophe, where parti- 
cles (stars) experience gravitational two-body interactions. 
The authors integrated the Boltzmann equation with a col- 
lision term due to gravitational two-body interactions, to 
obtain a set of moment equations. They then truncated, 
under reasonable assumptions, the hierarchy of such mo- 
ment equations such that only fluid-like conservation equa- 
tions remain. The effect of gravitational two-body interac- 
tions was naturally approximated as an effective heat con- 
ductio n. I n the CDM literature, iTevssicr. Chi eze fc Alirnil 
il997fl and lSubramanian. Cen fc Ostriker. 1200Q'1 followed a 
similar approach for the study of CDM haloes. Contrary 
to systems described by stellar dynamics where the num- 
ber of particles is small, gravitational two-body interactions 
are completely negligible for CDM haloes. The authors in- 
tegrated collisionless Boltzmann equation to obtain a set of 
moment equations, and truncated its hierarchy as done for 
stellar dynamics. This results in fluid conservation equations 
for collisionless systems. This idea may bother some readers 
since, strictly speaking, the collisionless nature of CDM pro- 
hibits the use of such an approximation. Collisionless parti- 
cles have, in principle, an infinite set of moment equations 
when the Boltzmann equation is integrated (BBGKY hier- 
archy; e.g. Birmcv &^Ireni<yjic 1982). However, a couple of 
simple assumptions enable us to treat CDM halo dynamics 
with the usual fiuid conservation equations. This section is 
devoted to the derivation of the fluid approximation, which 
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we will use in f^jto obtain self-similar equations both for the 
coUisionless CDM and the coUisional SIDM halo formation 
and evolution. 

I Alvarez. Ahn fc Shapirol ll2003ll have already applied 
this formalism to study CDM halo formation. Adopting the 
universal mass growth history r eported for CDM N-body 
haloes (e.g. IWechsler et al.ir2002h and applying the fluid ap- 
proximation, they have shown that the resulting shape of 
the equilibrium halo profile and its evolution match those 
of CDM N-body simulations remarkably well. In addition to 
this convincing result, we shall also justify the fluid approx- 
imation in iJ.S.ll 

We develop this model under certain conditions. First, 
spherical symmetry is assumed. The initial condition is given 
by a spherically symmetric overdense region, and the subse- 
quent evolution does not break the symmetry. Second, the 
infall of matter is assumed to be continuous. This can be 
achieved if we assume a smooth initial overdensity proflle. 
Third, we restrict our attention to the matter-dominated era 
such that the cosmic mean density pb oc t~^. This condition 
is true in the ACDM universe if the redshift is restricted 
to be l<z;<Zoq, where Zeq is the redshift for the matter- 
radiation equality. 



3.1 Fluid approximation of coUisionless CDM 
dynamics 

Let us describe the fluid approximation for a self-gravitating, 
weakly coUisional system in spherical symmetry. We define 
the average physical quantities as follows: 



d d d d 

Throughout this paper, we will use the Newtonian approxi- 
mation: when a system is much larger than its Schwarzschild 
radius and much smaller than the horizon size, motions of 
particles and the temperature of the system become non- 
relativistic. In this limit, we can use the non-relativistic 
Boltzmann transport equation. This equation can then be 
written more explicitly in spherical coordinates. For a sys- 
tem in spherical symmetry, / = /(|r|, v), and equation 1131 
reads 







dt ^ dr 



+ - (vl cot 9 - VrVe) 

r ^ ' OVg 

r dvd, 



dr I d 



(15) 



where $ satisfies the P oisson equation V^"I> — AnGp 
iBinnev fc Tremainell98'^ . By multiplying equation 1151 by 
J d^v v^Vg, where m, n are integer numbers, we can form a 
set of moment equations. Moment equations from the lowest 
order are 



dp ^ d 
dt r^dr 



(r^ipu)) = 0, 



+ ^(Pr + P'^^) + -{Pt - pe + pu^) 
dt dr r 

D ( pr\ du 

^:Dt U)+^'■aF^^^• 



(l6) 



= -p— , 17 
r^ 



(18) 





(7) 


j fcPy p J 


(8) 




(9) 


p{{Vr - {Vr)f), 


(10) 


p{{vg - {vg)f) = p{vl), 


(11) 


p{{v4, - {V4,)f) = p{vl), 


(12) 



D 



{A) 



Pr 



where / is the distribution function defined such that 
/(r, \)d^rd^v = mass within inflnitesimal volume d^rdPv 
at (r, v), p is the density, {A) is the average value of a cer- 
tain quantity A, u is the radial bulk velocity, Pr is the "ef- 
fective radial pressure", and pg is the "effective tangential 
pressure". Note that (yg) = {v^) = and pg — ptf, because 
of spherical symmetry. Anisotropy in the velocity dispersion 
occurs in general - i.e. pr 7^ pg or anisotropy parameter 
/3 7^ where f3 = 1 ~ y- ~ implying that we should treat 
Pr and pg separately. In a highly coUisional system, which is 
well described by fluid conservation equations, pr — pg and 
the usual pressure p = pr = pg. 

A self-gravitating system of coUisionless particles can 
be described by the coUisionless Boltzmann equation 

where ^ is the phase-space Lagrangian time-derivative, 
given by 



Dt \ 2p 



pg 



pgu 
r 



(19) 



where m is the mass enclosed by a shell at radius r, = 



at ^ "S' 



and 



Fl = ^{2{Vr-{Vr))vj) 



( 2 

2r2 dr 



'p{{V: 



Vr))vj)) . 



(20) 



(21) 



4r* dr 

Equations 11611 - 1191 1 are conservation equations of mass, 
momentum, "radial" energy, and angular momentum respec- 
tively. Note that equations II6II - 12111 are all in exact form, 
and the hierarchy of equations is not closed in principle. 

Now we make a further simplification that the distri- 
bution of Vr is skewless - vg and f0 are naturally skew- 
less because of spherical symmetry. In other words, we as- 
sume that Vr has a symmetric distribution around (vr). It 
is not straightforward to show that Fi and F2 are negligi- 
ble in equations 1181 and 1191 . However, we demonstrate 
that the assumption of "skewlessness" in the fiuid approxi- 
mation yields results which are in good agreement with the 
purely coUisionless CDM structure, for specific examples^. 



^ ISubramanian. Cen fc Ostrikerl l200(l) argue that an initially 
skewless distribution function / will remain skewless even after 
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-1 -0.5 
log,o (radius) 

Figure 4. Self-similar coUisionless halo fo rmation for e = 1: 
Phase-space diagram of iBertschingerl ll985^ solution as an illus- 
tration of the "skewless distribution" assumption. Along the line 
A there are a large number of shell-crossings, which are almost 
symmetric around (nr), and along the line C, as there is only one 
stream line, it is intrinsically skewless. On the contrary, there are 
only three stream lines along the line B, so the assumption is 
not good in the region inside and close to the outermost caustics. 
Tangential velocity distribution would have a similar behaviour 
when tangential motion is allowed. 



As shown in the IBertschingerl il985ft solution, for instance, 
coUisionless particles (shells) form a quasi-symmetric wind- 
ing structure in the phase space as shown in Figure 0] This 
indicates that we may neglect terms which arise as a result 
of skewness. Equations Hl()|l - l)18|l with the condition pg = Q 
and Ti = r2 = 0, 



| + 4;:('-^(H) = o, 



pu 



+ pu 



(22) 



Gm 



D 



Dt \2p 



Pr 



+ Pr 



du 
dr 



0, 



(24) 



can be used to solve purely radial problems, such as 
the spherical infall problems with similarity solutions by 
iBertschinge^ Jl985l e = 1) and lFiUmore fc Goldreichl Jl984l 
with e — 1/6). As seen in Figure^ the solution to equations 
12211 - 12411 is in good agreement with the true solution. We 
have also found, as seen in Figure |51 an excellent agreement 
between the coUisionless solution and the one produced by 
a radial fluid approximation in the case of e — 1/6 (the 



evolution. It is true if one is interested in the fine-grained distribu- 
tion, or each stream line. However, if one considers coarse-grained 
quantities as defined by equation JSJ, the contribution from mul- 
tiple stream lines should be counted and the overall distribution 
is not skewless, in general. 



£ = 1/6 case is of main interest in this paper). The differ- 
ence observed at caustics - places where the density becomes 
infinite - is negligible, because caustics do not affect the 
overall dynamics of the halo. Since the skew- free assump- 
tion naturally neglects dynamically unimportant structure 
(e.g. caustics) while accurately reproducing the profile of the 
exact solution in these radial cases, it may also be applied to 
describe CDM haloes, in which particles have a tangential 
motion as well. 

The final assumption is that inside the virialized 
structure (i.e. postshock region) the velocity dispersion is 
isotropic, or pr = pg. This is an empirical assumption: 
CDM haloes in cosmologi cal N-body simulation s show mild 
anisotropy. For instance, ICarlberg et al.l (Il997ll show that 
CDM haloes in their numerical simulation can be well-fitted 
by a fitting formula 



0{r) = f3„ 



4r 



+ 4' 



(25) 



where r is in units o f ^2 00, and Prn = [0.3 — 0.5] (see also, e.g., 
iThomas et al.ll998l and lColin. Klvoin fc Kravtsovl200nh . As 

we will show in the following sections, our similarity solu- 
tions have a virial radius r564 ~ 0.6r2oo, which then results 
in the maximum anisotropy f3 ~ [0.17 — 0.28]. This fact en- 
ables us to use /3 = to a good approximation. 

With these assumptions (spherical symmetry, skew-free 
velocity distribution, and isotropic velocity dispersion), the 
usual fluid conservation equations are obtained. They are 

d 



dt r^dr 



{r\pu)) = 0, 



d , s d , 2n 2 2 



p— 



DV2p' 



d 



■{r u), 



(26) 



(27) 



(28) 



p r^dr 

which are identical to the fluid conservation equations for 
a 7 = 5/3 gas in spherical symmetry. Resemblance of these 
equations to fluid equations indicate that we can expect an 
effective "shock" even for a coUisionless system because of 
the hyperbolicity of these equations. This is also illustrated 
in Figure 0] For instance, the density jump occurs when one 
moves from the "pre-shock" region (line C; one stream line) 
to the "post-shock" region (line B; three stream lines). How- 
ever, one should be careful in using this formalism because 
the approximation becomes worse where there are only a 
small number of phase-space windings. 

3.2 Fluid approximation of SIDM haloes 

We now consider the effect of SIDM coUisionality. We adopt 
a simple, heuristic approach to account for the heat conduc- 
tion as a result of finite cross-section. This will change the 
energy conservation equation (equ. 1281 ') to 



DV2p' 



d 

.2Qj. 



(r u) 



Vf, 



(29) 



where f is the heat fiux. 

The heat fiux resulting from SIDM interaction has been 
derived by BSI. We briefiy describe its derivation. When the 
scattering mean free path is much smaller than the system 
size (short mean free path limit, or diffusion limit), the heat 
flux (equivalent to in BSI) reads 
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■l-smfp 




(30) 



where b is an effective impact parameter of order unity, 
which is 1.002 for elastic scattering of hard spheres. On the 
contrary, when the mean free path is much larger than the 
system size (long mean free path limit), the proper length 
scale of heat transfer is not the mean free path (a multiple 
of the system size) but the system size, while the time scale 
is still the relaxation time. This results in 



Ilmfp 



3 , [p P d (p 

r^'^y-pi^^d-r [-p 



(31) 



where a = 2.26 for elastic scattering of hard spheres. 

In the end, we use a hybrid expression which is roughly 
valid in the intermediate regime as well as in two extreme 
regimes (short mean free path limit and long mean free path 
limit)^ 



-aba 




2 , 47rG 
aa H 



d_ 

dr 



(32) 



Throughout this paper, we adopt a = 2.26 and b — 1.002. 



3.3 Shock jump conditions 

As we have shown in previous sections, fluid conservation 
equations can be used and we expect an accretion "shock" 
to occur. Across the shock, the matter, momentum and en- 
ergy fluxes should all be continuous. Mathematically, a flux 
is a term acted upon by ^ in the conservation equations 
expressed in Eulerian coordinates. 

In the adiabatic case^, therefore, we obtain the usual 
adiabatic shock jump conditions from equations H26|l - H28|l : 



[pu] = 0, 

[p + pu'^] = 0, 



pu 



3p , 1-2 , P 

— It H — 

2p 2 p 



= 0, 



(33) 
(34) 

(35) 



where [A] = y4(preshock) — y4(postshock), u (= u — Us) is the 
bulk radial velocity in the shock frame, and equation 13511 
comes from equation 1281 if expressed in Eulerian coordi- 
nates by converting into ^ + 

For an SIDM the heat conduction is included, 

the energy jump condition will instead be 



3p ,1-2 , P , . 
2:^ + 2" +-p+-f 



pu 

where f — ff 



= 0, 



(36) 



^ Note that the intermediate regime described by equation 1321 
is indeed a mere interpolation of two different regimes. We be- 
lieve, however, that this is a good approximation: the intermediate 
value should not be too different from this smooth and continuous 
interpolation of two regimes. 

^ We will use the term "adiabatic" for the case where there is no 
heating mechanism (e.g. conductive heating due to elastic scat- 
tering) other than the shock heating. 



4 SELF-SIMILAR MODEL FOR SIDM HALOES 
IN THE MATTER-DOMINATED ERA 

In this section we flrst show that formation and evolution of 
SIDM haloes in galactic scale are well approximated by self- 
similar equations. We then apply the fluid approximation to 
the problem and derive a set of ordinary differential equa- 
tions. We will explain in detail how to solve these equations 
with proper boundary conditions. 

4.1 Self-similarity of SIDM haloes 

In H2.3I we described how a self-similar collapse model 
fits in cosmology by relating the logarithmic slope of the 
power spectrum, n, to the initial linear overdensity profile 
parametrized by e. We, therefore, should flrst know what n 
is relevant to the problem we solve. 

When coUisionality of particles enters the system, its 
corresponding length scale comes into play, which in gen- 
eral does not grow in proportion to rta. Similarly, a new 
time scale will enter the system. However, we can still make 
the system self-similar after the addition of coUisionality by 
some fine-tuning. The conductive heating term V • f enters 
the energy equation, and it should have the same time de- 
pendence as does the adiabatic change of the thermal energy, 
Pft (i)- Noting that Qc and V • f cx rL^-^ 

the condition 



C = 2, e 



(37) 



preserves self-similarity of the system. This condition, e = 
1/6, is equivalent to n = —2.5, which can be seen from 
equation Q. The turnaround radius and mass grow in time 
as 



(38) 



Now we find an intriguing coincidence. As described 
in i |2.3l n — —2.5 is a good approximate value for haloes 
of galactic mass. SIDM haloes in galactic scales are, there- 
fore, well described by self-similar equations. The following 
sections are dedicated to the detailed description of SIDM 
similarity solutions. 



4.2 Basic equations and problem solving scheme 

Under the condition of self-similarity, one can convert seem- 
ingly time-dependent equations into ordinary differential 
equations by properly scaling physical parameters. The 
turnaround radius rta is a natural choice for a length scale. 
With time t and the turnaround radius rta, we define di- 
mensionless physical quantities - radius, velocity, density, 
pressure, mass and heat fiux - as follows: 



A = r/rta, 

l/(A) = ^/(^), 

DiX) = p/pt, 



P{X)=P/ 



Pb 



M{\) = m/ {^^Pbrl 



(39) 
(40) 
(41) 
(42) 

(43) 



© 0000 RAS, MNRAS 000, 000-000 



10 K. Ahn & P. R. Shapiro 



FW = // 



Pb 



(44) 



Collapsing shells, which are assumed to be cold initially, 
obey Newton's law 



Gm 



(45) 



where m is the mass enclosed by radius r. For the initial 
density perturbation defined by equation Q, this Newto- 
nian motion can be describ ed by a se t of parametric equa- 
tions, as follows (Abadi, Bower fc Navarro. .20001 see also 
lBertschingeilll985l for case of £ = 1): 



A : 



\e/2) 



V(\) = A 
D{\) = 



(1 - cos 61)2 
9 {6 - sin ( 



2(1 



s6')3(l + 3ex)' 



(46) 
(47) 

(48) 

2(l-cose)3' ^^^^ 

where x = l-(3/2)(K(A)/A). Equations ^ - ^ therefore 
describe preshock motion. 

The postshock motion of shells is described by full hy- 
drodynamic equations. Equations I26II . 1271 . 1291 . 1321 and 
the definition of the infinitesimal mass dm — 4-Kr^dr can 
be written with these dimensionless quantities as a set of 
ordinary differential equations: 

(V - 2A)D' + DV' + - 2D = 0, 



A 



{V - 2X)V' + V = 



EL 

D 



2 M 
9 A2' 



(F-2A)(^-^U-H 



2(FA')' 
3PA2 



F = 



Zah 2 




M' = 3A^_D, 



(50) 
(51) 
(52) 

(53) 
(54) 



where the prime indicates differentiation with respect to A, 
and the nondimensional coUisionality parameter Q' is de- 
fined as Q' = CTpfcT-ta. We will later use a nondimensional 
constant Q = ap^rs — X^Q' , which is more directly related 
to the collision rate of an SIDM particle in a virialized struc- 
ture with the shock radius r^. For instance, the number of 
collisions a particl e experiences in a time At is given by 
N = apVreiAt fe.g. lBurker3l2000l) . The conversion of Q into 
A'^ is straightforward: 



N 



p Urol At 



p p At 



Pb 



P Pb Ts 



(55) 



Here we used the relation v^ci = o,yP/Pj where a = 2.26 
again, which relates the average thermal velocity to the 
relative velocity for particles interacting elastically as hard 
spheres (equs. [7.10.13], [12.2.12] in iReif 1965.1 . 

Different solutions arise for different values of Q{Q'). 
To solve the coupled ordinary differential equations (equs. 
HO] - ED)) we need to connect the preshock values given 



by equations 1461 - 1491 with ^ = 2) to the postshock val- 
ues. This is described by the shock jump conditions (equs. 
1331 . 1341 . and 1361 ') and the continuity of mass, which are 
expressed by nondimensional variables as 



D2 = 



2F2 



P2{Vi-2K 



+ 



2P2(Vi -2As) 



4Di, 



D^{Vi-2K 



V2 = 2K 



Ma = Ml, 



2F2 



P2(Vi -2A, 



2A,), 



(56) 



(57) 



(58) 
(59) 



where the subscript 1 denotes preshock values, while 2 de- 
notes postshock values. Note that Pi = and Fi = be- 
cause we assume cold infall. Note also that without terms 
containing F2, equations 1561 - 1591 are identical to the adi- 
ab atic jump conditions for ^ = 5/ 3 gas (see equs. [6a] - [6d] 
in I Abadi. Bower fc Navarro! I2OO0I) . These additional terms 
arise because of finite conductivity in the postshock region. 
Finally, the inner boundary conditions are 

M{\ = 0) = 0, (60) 

V(A = 0)=0, (61) 

and 

F(A = 0) = 0. (62) 

In principle, the fluid equations (equs. 1501 - 1541 '). 
preshock equations (equs. 1461 - \49\ ). jump conditions (equs. 
1561 - \b9\ ). and inner boundary conditions (equs. 1601 - 1621 ') 

yield a unique solution. This solution was obtained numeri- 
cally, by iteration, as follows. We arbitrarily chose the shock 
location As (close to that of the adiabatic solution), central 
density D{0) and central pressure P(0). We then integrated 
differential equations from A = outward, using the LSODE 
(Livermore Solver for Ordinary Differential Equations), in 
which case variables behaved well. At the chosen shock lo- 
cation, we then obtained preshock values using the jump 
condition. If these values differed from those obtained from 
equations 1461 - 1491 . we went back and chose a different 
As, -D(O), and P(0). We iterated this process until the jump 
condition was satisfied to a given error tolerance. This way 
we could get approximate solutions for several selected val- 
ues of Q (listed in Table 0, all of which satisfied the shock 
jump conditions (equs. 1561 - 1591 ') to less than 0.1% error. 
This was a very tedious and time-consuming job, but au- 
tomation of the iterations made it possible to achieve this 
goal. 

In practice we could not perform integration from A = 
because the differential equations have a coordinate singu- 
larity at the centre. Instead, we performed the integration 
from some small A, using asymptotic forms for A ^ 1: 



D 

F : 



D{0), P 
-|P(0)A. 



P(0), V 



A, M 



D(0)A^ and 



5 RESULTS 

We solved the coupled set of ordinary differential equations 
which yield the similarity solutions for SIDM haloes for the 
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Figure 5. Similarity solution dimensionless profiles for low-Q 
regime. Q = means "no conduction," i.e. "adiabatic" post-shock 
gas. As Q increases, core density decreases and core temperature 
in creases. Dimensionle ss similarity variables follow the definitions 
in lBertschingeJ (^^. 



Figure 6. Similarity solution dimensionless profiles for the high- 
Q regime. Profiles are indistinguishable from those in Figure 151 
even though the Q values are quite different. The effect of Q 
is reversed compared to the low-Q regime: as Q increases, core 
density increases and core temperature decreases. 



full range of values of the dimensionless collisionality param- 
eter, Q. The results are plotted in Figures |K| and |S| Before 
we discuss the solutions for SIDM haloes with conduction 
{Q 7^ 0), we begin with a description of the solution for 
(5 = 0, the case without conduction. As we shall see, the 
latter yields a density profile which is in good agreement 
with the results of N-body simulations of haloes in the stan- 
dard CDM model. We are justified, therefore, in applying 
the same self-similar infall model with Q 7^ to describe 
SIDM haloes. 

Let us briefly describe the general properties of the 
SIDM similarity solutions and emphasize the importance of 
cosmological infall, before we describe the results in full de- 
tail. 



5.1 Self-similar haloes without 

conduction(e — 1/6): an analytical model for 
CDM N-body results 

Before treating soft-core solutions, we describe properties 
of the adiabatic solution for e = 1/6. As we shall show in 
what follows, the adiabatic infall solution for e = 1/6 case 
resembles standard CDM haloes in many respects. First, it 
has a density cusp with a logarithmic slope ~ —1.27 for 
4 X 10~^ < r/r2oo < 1.4 x 10~^, where r2oo is the radius in 
which the average density is 200pf,, if the density is extrap- 
olated beyond Va with the NFW profile that best-fits this 
adiabatic solution (Fig. 0. Note that it is possible to have 
a slope shallower than —2 because particle velocities are al- 
lowed to have a tangential component. If no tangential mo- 
tion is allowed, the value —2 is the s hallowest slope possible 
for coUisionless dark matter haloes jRichstone & Tremainel 



m. 

19841: iTevssier. Chieze fc Alimilll997l - TBcrtschingcr 1998) 



Second, the temperature profile is very similar to that 
of CDM haloes. The temperature is zero at the centre, rises 
to a maximum at some A, and then falls to a nonzero value 
at the shock. The most important part is that the tempera- 
ture inversion exists. If not, the addition of conductivity will 
worsen the situation by initiating g ravothermal catas trophe, 
rather than generating a soft core iColm et al.ll20o3) . 

The average density inside the shock radius is found to 
be 564pi, for the adiabatic solution. This value is larger than 
the average density (~ 200) usually adopted by convention 
to identify virialized CDM haloes in N-body simulations. In 
terms of radius, Vs = r564 — 0.6r2oo. Our similarity solu- 
tion, therefore, yields a shock at a smaller value of radius 
than is typically used to characterize the virial radius of 
CDM haloes in N-body simulations. However, as long as 
we focus on the properties inside r564, the adiabatic solu- 
tion is a good fit to CDM haloes as we shall see below. 
As seen in Figure |7| except for the innermost central re- 
gion whCTe_th«_ambi£i^^ slope e xists between 
-1 llNavarro. Frenk fc WhitellT997i) and -1.5 iMoore et alJ 
Il999h . our adiabatic solution agrees with the NFW profile 
and the Moore profile to within 10%, depending on the 
concentration parameter. Compared to halo profiles stud- 
ied by Dicrnand, Moor e fc Stadel (2Q04), the agreement is 
even better (Fig. |7|). We also find that a local logarith- 
mic slope slowly changes to shallower values as one ap- 
pr oaches the centre, wh ich agrees with the trend reported 
bv lNavarro et all ll2004ll . We demonstrate this as follows. 

We compare the adiabatic solution with the CDM N- 
body halo profiles mentioned above: i.e. the NFW profile, 

Pnfw = — 2 ' (^^) 

(r/rsc) (1 + r/rac) 

the Moore profile. 
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Figure 7. Comparison of self-similar halo profile without con- 
duction with N-body results for CDM haloes: (top panel) Den- 
sity (in units of cosmic mean density) vs. radius (in units of cur- 
rent turnaround radius) for similarity solution (Q = 0; e = 1/6) 
(solid), the best-fitting NFW profile (c = 3), Moore profile 
(c = 1.6) and several a/37 profiles (« = 1, /3 = 3) for dif- 
ferent values of 7, 7 = 1 . 3 (be st fit), 1.02, 1.16, as labelled. 
iDiemand. Moore fc Stadell <2004h find N-body results for CDM 
haloes best fit by (a, /3, 7) = (1, 3, 7 = 1.16 ± 0.14) (shaded re- 
gion); (bottom panel) fractional deviation of the N-body results 
from similarity solution. Note that A = As ~ 0.09 corresponds to 
J-s = r-564 ~ 0.6r200- 



(64) 



(r/rec)^'^ (1 + (r/r.c)^-') ' 
and the a/37 profile (e.g. iDiemand. Moore fc StadellEo04l) 



Pa(3^ — 



{r/r^^y (1 + (r/rsc)") 



(/3-7)/c« ' 



(65) 



by finding the best-fitting parameters. The NFW and Moore 
profiles have two free parameters (psc and r^^), while the 
Q/37 profi le has three (psc, t'sc, and 7) wh en a and /3 are 
fixed as in lPiemand. Moore fc Stadell (|2004l ). We depict this 
comparison in Figure|7| with the "concentration parameter" 
- defined by c = r2oo/''sc - used to find the best-fitting N- 
body halo profiles. For q/37 profiles, wo use the best-fitting 
profiles by Dicmand, Moore & Stadcl (200^, namely a — 1, 
/3 = 3, 7 = 1.16±0.14. Among these profiles, the a/37 profile 
with a = 1, f3 = ?> and 7 = 1.3 provides the best agreement 
with the adiabatic solution. 

The relatively low concentration parameter, c ~ 3, re- 
quired for the best-fitting NFW profile deserves attention 
with respect to the cosmological mass accretion rate. This 
value is a bit small compared to the typical concentra- 
tion parameters, c ~ [4 — 20], observed at 2: = in N- 
body simulations. Recent high resolution N-body simula- 
tion results, however, report such low concentration param- 
eters if CDM haloes a re observed at higher redshifts (e.g. 
iTasitsiomi et alJl2004l) . Individual haloes in CDM N-body 
simulations evolve over time, on average, through a contin- 



Figure 8. Profiles of heat flux and bulk velocity. Negative heat 
flux in the core indicates that heat is transferred from the halo 
to the core. The positive velocity bump indicates that the central 
density keeps being flattened by a net expansion of the core. Two 
different Q values were selected for comparison. 



uous sequence of universal-shaped mass profiles of increas- 
ing total mass ([ Tasitsio mi et al.ll2004t Ivan den Boschir2002l : 
IWechsler et al.ll2002^ . This Lagrangian mass evolution can 
be characterized by a universal mass accretion history: 



M{a) = Moo exp(-2a//a). 



(66) 



where a is the cosmic scale factor and aj is some partic- 
ular value of a, such as that at which dlogAf/dloga = 2 
iWechsler et al.l l2002l) . As the mass of each halo grows 
with time due to the average effect of mergers and smooth 
infall, so does the concentration parameter c of its den- 
sity profile, roughly as c{a)/c{af) oc a/a^ for a/a{ > 1 
JWechsler et al.l2002h . after hovering at low values c ~ 2 — 4 
during the initial phase of most rapid mass assembly prior 
to Uf llTasitsiomi et all [2004'). Our similarity solution has 
''diog^a ~ ^' seen in equation 1381 . and this corresponds 
to a = a//3. As our solution corresponds to a very early 
epoch in the halo formation history, or a fast accretion rate, 
such a low concentration parameter is a natural outcome. 

In summary, we have shown that the adiabatic solution 
for £ = 1/6 approximates the N-body CDM halo profiles 
well, thus providing physical insights about their origin. In 
the foUowing sections, we will describe SIDM halo solutions 
which arise as a result of non-zero conductivity. 



5.2 Self-similar haloes with conduction(e — 1/6): 
an analytical Model for SIDM haloes 

5.2.1 Low-Q regime 

We define the low-Q regime as Q < Qth ~ 7.35 x 10~*. 
All the solutions have an isothermal, flat-density core, ex- 
cept for extremely small Q, where the system undergoes 



© 0000 RAS, MNRAS 000, 000-000 



Self-interacting dark matter haloes 13 



Table 1. Central parameters and the shock location for different 
Q solutions. This table lists the whole range of Q. 
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1.062{- 


■1) 


9.133(- 


-2) 


oo 




oo 


oo 







9.0434(- 


-2) 



adiabatic collapse just as it would in the absence of SIDM 
conductivity. In this regime, as Q increases, the core density 
and pressure decrease. In other words, higher SIDM coUi- 
sionality corresponds to a flatter core. This trend is also 
observed in the temperature profile. As Q increases, heat is 
more effectively transferred into the centre to equalize the 
temperature (Fig. Dependence of the central density, 
temperature and the location of the shock on Q is listed in 
Table [U (see also Fig. 0). 

This regime roughly corresponds to the long mean free 
path limit. In this case, the heat fiux is an increasing function 
of the SIDM cross section a (eq. 1311 ). As an increase in Q 
is achieved by an increase in cr, the heat flux also increases 
correspondingly. In a more detailed description will be 
given about the quantitative relation between the low(high)- 
Q regime to the long(short) mean free path limit. 

In H5.5.2I we will also show that cosmological SIDM 
N-body simulations to date have been performed only in 
the low-Q regime. They see a monotonic behaviour of the 
halo profile depending on cr - as cr increases, core density 
decreases. In the next section, we will show that there is an 
additional regime, the high-Q regime, where this behaviour 
is reversed. 

5.2.2 High-Q regime 

In the high Q regime, solutions have Q > Qth = 7.35 x 10"'*. 
Once again, all the solutions have an isothermal, flat-density 
core, except for an extremely high Q case. In the high-Q 
regime, however, as Q (or cr) increases, the core density in- 
creases (Fig. Table 0, contrary to the behaviour ob- 
served for the low-Q regime. Therefore, Qth gives the solu- 
tion with the minimum possible core density pcoro — W^pb- 

This behaviour occurs because the core of the halo is 
now in the short mean free path regime. In this case, con- 
trary to the low-Q regime, the hybrid conduction term, equa- 
tion I32II . converges to the expression valid in the short mean 
free path limit, equation 13UII . An increase in Q or cr, there- 
fore, results in a decrease in f. Physically, the mean free path 
decreases as Q or cr increase, and it results in reducing the 
heat conduction. 

For an extremely high Q (or cr), we found that the 
solution becomes identical to that of the adiabatic infall 
case. This result agr ees qualitatively with the result by 
lYoshida et"ai] ll2000al) and lMoore et alj ^00), where they 
performed a smoothed particle hydrodynamics simulation, 
corresponding to the fluid limit of an infinite cross section. 



What they found in haloes was a density cusp instead of a 
flat-density core. Quantitatively, our result does not fully 
agree with their result in which they find a profile even 
steeper than that of the coUisionless case. As pointed out 
bv iYoshida et al . (2000aJ , this may be attributed to the fact 
that small-scale shocks arise in the case of the fluid regime, 
thus increasing the entropy and ultimately steepening the 
central density. Our model is based upon an assumption that 
mass accretion is smooth and, therefore, cannot reproduce 
this effect. 



5.3 Meaning of the collisionality parameter Q 

We have showed that solutions are parametrized by the colli- 
sionality parameter Q, and also that there exist two regimes 
divided by a threshold value Qth- To understand the physi- 
cal meaning of Q, we now describe two relevant quantities: 
mean free path and the number of scatterings an SIDM par- 
ticle experiences per Hubble time. 

We first show that the ratio of the mean free path to 
the gravitational scale height (to be explained below) in the 
centre is very closely related to Q, and it provides a very 
clean explanation of the behaviour of the similarity solu- 
tion. In previous sections, we explained the opposite trends 
observed in the two different regimes in terms of the heat 
flux f (equ. 1321 ). The dependence of f is determined by the 
ratio 1]^ = {4TrG/p)/{aa^): f oc cr for i]^ ^ 1, and f oc 1/a 
for ri^ <^ 1. We find that, indeed, ri is the ratio of two length 
scales, the mean free path and the gravitational scale height. 
The gravitational scale height (BSI), 

H = ^a^/i4nGp) = Vp/(4^Gp2), (67) 

is a rough measure of the size of a given self-gravitating sys- 
tem, where ay is the velocity dispersion - also note, how- 
ever, that equation 16711 determines the local gravitational 
scale height. In terms of the mean free path Amfp = l/(pcr) 
and a = 2.26 as in equation (I31II . ri = Xmip/iVaH). We 
find that t; = 1 at the centre for Q — Qth. rj at the centre 
monotonically increases as Q increases. The low-Q regime 
then corresponds to a condition rj > 1, and the high-Q to a 
condition r] <1. Dependence of 77 on Q, as well as its radial 
variance, is plotted in Figure |5] 

The number of scatterings that an SIDM particle expe- 
riences during the age of the universe, which we denote by 
A'^, is also an interesting quantity. As higher Q means more 
frequent scattering, A'^ also shows a monotonic dependence 
on Q (Fig. 1101 ) as rj does. According to equation I55II . it 
can be shown that N is expressible in terms of dimension- 
less quantities as 

where the subscript "0,th" refers to the value at r = for 
Q = Qth, and Ao^th = 129. It is interesting to see that 
Ao.th ~ 100 is required to achieve the maximal conductivity, 
namely Q = Qth- We will handle the significance of this 
quantity in >)5.5.2I when we compare our result to SIDM 
N-body simulation results. 

Finally, we stress the importance of the radial variance 
of ry and N . When Q/Qth <C 1 or Q /Qth S> 1, the system can 
be said, in a global sense, to reside in the long mean free path 
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Figure 9. The parameter q versus dimensionless radius A, where 
T) = ^^{p/{y/aH), the ratio of mean free path to gravitational 
scale height, for different values of Q. From top to bottom, each 
curve corresponds to Q/Qth = 1-23 • IQ-^, 9.77 ■ 10"^, 6.16 • 
IQ-^, 1, 17, 59.4, and 186, respectively. The constant a = 2.26. 




Figure 10. Number of scattering that an SIDM particle expe- 
riences during the age of the universe in dimensionless radius 
A. From bottom to top, each curve corresponds to Q/Qth = 
1.23 ■ 10^3, 9.77 ■ 10-3, 6.16 • IQ-^, 1, 17, 59.4, and 186, respec- 
tively. 



limit or in the short mean free path limit, respectively. When 
Q/Qth ~ 1, however, such a global definition is not valid. In 
the Q = Qth solution, for instance, 77 = 1 at the centre while 
?7 « 10 at the shock. A'' varies from ~ 100 at the centre to ~ 2 
at the shock. This example shows that a global assumption 
of the short (or long) mean free path limit is not always 
valid, which requires a more careful attention when Q ~ Qth. 
We will handle this issue again in ^ with respect to the 
estimate on the ram-pressur e stripping of substructure in a 
cluster environment ma de bv lFurlanetto fc Loebl i2002f) and 
iNataraian et alJ i200'i) . 



5.4 Importance of cosmological infall 

The SIDM core grows in size as a fixed fraction of the 
turnaround radius, as guaranteed from the beginning be- 
cause of the self-similarity of the system. 

Therefore, this model shows that cosmological infall at 
a certain rate can completely inhibit the gravothermal catas- 
trophe of the core by constantly pumping hot material into 
the halo (see Fig. This contradicts the prediction made 
in previous studies that the core would suffer g ravothermal 
catastrophe in a Hubble time fe.g. lBurkertl200 [ih. Even when 
the infall rate is smaller than that required for our similarity 
solution (Af oc f'*), it will inhibit the core collapse to some 
extent. Only when the infall rate drops to a point where a 
system can be considered isolated will previous estimates of 
the timescale of collapse be valid. The net effect is a sub- 
stantial delay of the collapse phase. 



5.5 Application 

5.5.1 Collisionality parameter as a function of a and M 

So far, we have described solutions in terms of Q. Now we 
seek a way to apply our solutions to practical problems. This 
is done by obtaining the dependence of Q on the scattering 
cross section a and the mass of haloes Ai. From the Press- 
Schechter formalism, we find Q = Q{(t, M) for "typical" 
haloes of mass M, which collapse when om = (5crit (cta/ is the 
standard deviation of the density fiuctuations at the collapse 
epoch 2coii(M) according to linear perturbation theory after 
the density field is filtered on the scale M; (5crit is the value 
of overdensity linearly extrapolated to a moment when the 
nonlinear overdensity becomes infinite.) Usually these are 
called l-o"M fluctuations. 

As the e = 1/6 adiabatic infall solution has a shock at 
r564, the mass contained inside 7-564, M564, is smaller than 
M200, which is typically quoted in the literature. In other 
words, the shock location is displaced from r2oo substan- 
tially. As ra is a function of M and Zcoii, in order to get 
Q = Q{M, Zcoii), we first should relate M564 to M. We there- 
fore need a model whose density profile extends at least to 
7'2oo . Wo use the "truncated isothermal sphere " (TI S) model 
( Shapiro, Iliov & Raga 1999; Ilicv & S hapirollioOlh for this 
purpose, for the following reasons: (1) it has a unique den- 
sity profile, (2) all physical quantities are fully determined 
by the values of M and Zcoii, (3) it has been proven to agree 
well with CDM prediction in many aspects IS hapiro fc IlievI 
1^002') ■ and (4) the average quantity (i.e. the average temper- 
ature) inside its own r^ei is in a good agreement with that 
of our similarity solutions and of CDM N-b ody haloes^ . We 
use a convenient set of formulae given by llliev fc Shapiro! 
i200ll) to get rs = r^etiM) for a given halo of mass M at 
its collapse epoch. In this model, the mass of the halo M is 
enclosed by the "truncation radius" rt, given by 



rt = 187.2 



M 



1012/i-iM, 







1/3 



!^//^(l+2coii)"'/i"'kpc, (69) 



* We find the same level of agreement between our similarity 
solutions and the N- body haloes as was found bet ween the TIS 
and N-body haloes in lShaoiro. Iliev fc R.agal il99HL see the mass- 
temperature relation in § 8.4). 
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and we find that 
M564 = 0.587M, 
and 

r564 = 0.514 rt 
= 96.22 



(70) 



(10^™^) "o^^^(l + ^con)-/.-kpc.(71) 

Note that M in this model is equivalent to M130 . Equations 
16911 ■ 17011 ■ and 17111 are valid only when the universe is in 
the matter-dominated era. The mean matter density is in 
general given by 

Pb{z) = f7oPo,crit(l+z)^ = 1.88xlO-'''f7o(l+^)^/i'g/cm^(72) 

When the scattering cross section a is a constant, Q also 
remains constant in a matter-dominated era where pt oc t^^ 
if £ = 1/6 (or Ts oc t'^). From the equations above, we find 
that 



Q 



(^) 
M /I 

0.7^ Vl 



^coll 



218.5 cm2g- 
2 



M 



(lOio/i-iMQ 



1/3 



2.09 



Note that typical {1-um fluctuation) haloes of M = 
10^°/i~^Mq collapse at Zcoii ~ 2.09 in the currently- 
favoured ACDM universe with h = 0.7, = 0.27, = 
0.73 and as = 0.9. If we restrict the mass range of haloes 
such that the high mass end will still collapse in the matter- 
dominated era - «coii ^ 1 ~ and the low mass end roughly sat- 
isfies the self-similarity, we should choose haloes with masses 
M ~ [10® - 10^^] h~^MQ. In this mass range, equation ff^ 
reads 



9.31] X 10" 



218.5 cm2g-i 



(74) 



for a ACDM universe (see Fig. for a = 218.5 cm^g"-^, 
Q = 1.68 X 10"* corresponds to M = lO^h'^MQ, while 
Q = 9.31 X 10"" corresponds to M = lO" ''^/i"^M0). Here 
we applied the Press- Schechter formalism to obtain Zcoii(Af). 
However, just for comparison, we also calculated Q for haloes 
with M> 10^^ Mq , which typically have collapsed only 
recently. Such massive haloes must be rare, high-cr fiuctua- 
tions in order to collapse by the present, and we cannot ap- 
ply our similarity model to these haloes because their mass 
assembly deviates substantially from self-similar evolution. 
Note that the rarer the objects are, the higher the collapse 
redshift 2coii is, which in turn makes Q larger in equation 
1731 . This trend can be seen in Figure [TTI 

5.5.2 Comparison with N-body simulation results 

Consider the range of a id entified in the cos mological N- 
body simulation of SIDM bv lDave et all ll200j) : a = [0.56- 
5.6] cm^g"^ . In a ACDM universe, this range of a in equation 
1741 then yields the range 

Q ~ [4.31 X 10"^ - 2.39 x 10"^] 

~ [5.86 X 10"* - 3.25 x 10"^] Qth (75) 

for haloes with W^H-^Mq < M < IO^^/i-^A^q. Here 
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Figure 11. Left; Q vs. mass of haloes at their typical formation 
epoch for different a. From bottom to top, curves correspond 
to CT = 0.56, 5.6, 218.5, 1.2 x 10*, 2.7 x 10* respectively. They 
all correspond to I-ctjv/ density peaks (ctm means the standard 
deviation of the density fluctuations filtered on mass scale M at 
the collapse epoch); Right: Q vs. mass of haloes at their formation 
epoch for v-crm (v=l, 2, 3) fluctuations with a = 218.5 cm^g~^. 
Cluster-sized haloes observed at present will be clustered around 
the crossing point of the S-ctm Hu^ Eind 2:t,oii = line. 



M 



IO^H-^Mq, while Q = 2.39 x 10"^ corresponds to 



(73) a = 5.6 cm^g-^ and M = lO^^_^Hr^J^ 



The simulation results o f iDave et al] ll200lfl reside in the 
low-Q regime, as seen in equation 1751 . Our solutions allow 
us to identify a corresponding range of Q- values in the high- 
Q regime for which the density profiles are indistinguishable 
from their low-Q counterparts. The range of solutions de- 
scribed by equation 1751 can be matched by solutions with 



[1.37 X 10"^ - 0.17] ~ [18.6 - 231] Qth, 



(76) 



or (J ~ [1.2 X 10* - 2.7 X 10*] cm^g'^ Here Q = 1.37 x 10"^ 
corresponds to cr = 1.2 x 10*cm^g"^ and M = IG^H'^Mq, 
while Q = 2.1 X 10"^ corresponds to cr = 2.7 x 10* cm^g"^ 
and M = IO^^/i'^Mq. 

Therefore, we predict that a cosmological N-body sim- 



produce results similar to those of the iDave et al.j (12001 
simulations. The nondimensional core density -D(O) corre- 
sponding to the range of Q found in this section lies in 
the range Z)(0) ~ [3.3 x 10* - 7.4 x 10^], while the nondi- 
mensional core temperature Tcorc = is in the range 
Tcoro — [0.11 — 0.41]. However, in order to obtain observa- 
tionally acceptable values, we should actually fit the implied 
rotation curves of our similarity solutions directly to the em- 
pirical data. This is the main topic of H5.5.3I 



5.5.3 Rotation curve fitting 

In this section, we find best-fitting similarity solutions which 
match the observed rotation curves of dwarfs and LSBs. To 
do so, we simply compare our similarity solutions to an em- 
pirical fit bv iBurkcrt, 11995.') . Burkert found that a density 
profile given by 



p{r) 



pod 



(r-t-r-o)(r2 -f rg)' 



(77) 



Q 



4.31 X 10 ^ corresponds to a = 0.56 cm^g ^ and 



where po and ro are free parameters which represent the 
central density and a scale radius, respectively, matches halo 
density profiles which are derived from the observed rotation 
curves of dwarf galaxies. Recent studies of high-resolution 
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Table 2. Best-fitting concentration parameter of the Burkert pro- 
file and /v. x^/v is normafized by the value found for Qth 
solution, (x'^/^)Qth ~ 2-06 X lO"**, the minimum. 



Q/Qth 


'^564, Burkert 


(xV^')/(xV^)q,, 





5.51 


4.72(2) 


1.23(-3) 


5.51 


3.86(2) 


9.77(-3) 


5.49 


1.36(2) 


6.16{-2) 


4.98 


2.51(1) 


1 


3.95 


1 


1.70(1) 


5.17 


1.11(1) 


5.95(1) 


5.64 


9.82(1) 


1.86(2) 


5.59 


2.99(2) 


oo 


5.51 


4.72(2) 



Hq rotation curves of d warfs and LSBs confirm this result 
iMarchesini et al .120021: see also references therein): they use 
a "hybrid" of Ha and HI rotation curves which can extend 
from the centre to rmax, and find that the Burkert profile is 
the best fit for this range. 

For fitting purposes, one should specify a radius where 
the circular velocity [v = Gm/r) is the same for different 
halo models. In other words, local density may vary but the 
mass enclosed by such a radius - let us denote it by the 
"normalization radius" Vn - should be the same. We choose 
fn = f564, the "shock radius" of our similarity solutions, to 
find the best-fitting similarity solution to the Burkert profile. 
In this case, the "concentration parameter" 0554, Burkert = 
r^m/'ro is the only free parameter. The goodness of a fit is 
observed through the relative mean square deviation, 

= V ( ^(''O- "Burkert (rO \ ' ^^^^ 

V "i^Burkortfri) / 

i ^ ^ 

where is the radius of the ith data point and is the 
number of such points. 

We find that the solution with Q — Qth is best fit to the 
Burkert profile, as seen clearly in Figure [Hand Table 12'^''. 
As shown in the right panel of Figure^] when Q — Qth, the 
SIDM halo rotation curve is virtually indistinguishable from 
the Burkert profile at all radii r < rmax, by contrast with 
the strong disagreement at small radii between the Burkert 
profile and the NFW and Moore profiles, respectively. Since 
the Q = Qth solution has the most effective conductivity, 
which is closest to an isothermal structure among the sim- 
ilarity solutions, we argue that dwarfs and LSBs described 
by the Burkert profile are systems which are almost fully 

^ We find the same answer, that the Qth solution is the best-fit 
to the Burkert profile, when we relax the constraint r„ = r5g4 
and assume that the two parameters, po and ro, are both free. 

The values of /'^ shown in Tablel^show the relative quality 
of the fits of the Burkert profile to our self-similar SIDM profiles 
for different values of Q. In order to interpret these values in an 
absolute sense to determine if the fits are acceptable, we need 
to know the uncertainty of the fit of the Burkert profile to the 
observational data points. For example, if the l-cr error bars are 
all a fraction / of the value of the data points, then the quantity 
/"■^(x^/^) should be small, i.e. < 1, for the theoretical curves to 
be a good fit at the l-cr level. Since (x^/!^)Qth = 2. 06 X 10~*, 
the SIDM profile for Q = Qth is a good fit to the data as long as 
0.014. 



relaxed. This argument is supported by the fact that the 
TIS solution is almos t identical to the Bur kert profile (see 
Fig. [12]; also refer to llliev fc Shapirdl200lD . The TIS solu- 
tion is obtained by assuming that the system has a uniform 
temperature, isotropic random velocity, and the minimum 
possible energy, which in effect is equivalent to assuming a 
fully relaxed system. The solution with Q = Qth corresponds 
to the most relaxed system among our similarity solutions. 

5.5.4 High value of a: Contradiction with SIDM N-body 
simulation results? 

An attempt to identify the range of SIDM cross section re- 
quired to produce density profiles in agreement with dwarf 
galaxy rotation curves was previously made using N-body 
simulations of SIDM halo formation from Gaussian-random- 
noise cosm ological density fiuctuations l|Dave et alJ I2001I : 
lYoshida "et al., 2000b) . These N-body results also indicated 
that, for this range of SIDM cross section, larger mass haloes 
(i.e. from the Milky Way to clusters) produce density profiles 
with flattened cores which are even more pronounced than 
those of dwarf galaxies. Since astronomical evidence suggests 
that such large haloes have relatively smaller cores, if any, 
than dwarf galaxies, this has led to the suggestion that the 
SIDM cross section must depend upon the relative velocity 
of the scattering events, decreasing with incre asing veloc- 
ity to suppress this effect in larger haloes (e.g. IColm et alJ 
l2002h . How do our self-similar solutions for SIDM halo for- 
mation compare with these N-body results? 

The value of Q in the similarity solutions which best 
fits the dwarf galaxy rotation curves, Qth, corresponds to 
a value of the SIDM cross section when we identify the 
halo mass and redshift to which we apply the similarity so- 
lution, as described above in § 15.5.11 For haloes of mass 
M ~ 10^'^ Mq, which roughly represents the mass of 
dwarfs and LSBs observed, the solution with Q = Qth im- 
plies that a ~ 218.5 cm^g~^, if the observed galaxies formed 
at the typical epoch for their mass scale (i.e. I-ctm fluctu- 
ations). This value is significantly larger than the range of 
acceptable cross section values reported for the N-body re- 
sults for SIDM h aloes bvlDave et alJ J200lf) (equ. [73) and 
similar results bv lYoshida et alJ J2000bl) ! 

Such discrepancy is observed also in A'^, the number of 
scatterings that an SIDM particl e exper iences during the 
age of the universe. 'Yo shida et all J2000bl) report that A*' ~ 
1— 10 is enough to generate a soft core in their N-body SIDM 
haloes. We have seen in § 15.5.31 however, that Q « Qth is 
required to find acceptably flattened soft cores, and this in 
turn corresponds to A*' ~ 100 in the core region as described 
in g lO 

We may attr ibute this discrepancy to the fact that 
iDave et alj i200 j) did not actua lly perform rotatio n curve 
fitting with their N-body results. iDave et al.l ll200ll) instead 
found their preferred value of a by constraining the halo 
density at r = Ikpc for haloes at the present epoch to be 
in the range [0.01 — O.IJMqpc"''. However, limited numer- 
ical resolution prevented them from determining the halo 
density at radii as small as those required to match the ob- 
served rotation curves which show fiat-density cores. Our 
results suggest that, had their simulations been capable 
of resolving the profile at smaller radii, they would have 
found that the density continued to rise to a higher value 
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Figure 12. Rotation curve fitting. The upper left panel compares several different Q-solutions to the Burkert profile, normalized to r^gn. 
The lower left panel shows the best-fitting solution to the Burkert profile, namely the Qtjj solution. The right panel compares various 
halo models to the Burkert profile, normalized to rmax. It has the Q^i^ profile in r and v in units of r„,ax Burkert = 0.835 r564 giDM and 
''max, Burkert = 1-01 f564,SlDMi respectively, for the same profile as plotted in the lower left panel. In both boxes, the top panel shows the 
relative difference of a given profile, [v — fBurkcrt ) /''Burkert ■ The line types of the upper right panel follows the meaning of those in the 
lower right panel. 



at smaller radii. As seen in Figure 1 of iDave et alJ i200ll) . 
three haloes with M ~ 10^ Mq (M = (1.7,0.9, 1.1)1O^M0) 
are almost unaffected by the inclusion of SIDM coUisionality 
ffcr ~ [0.56-5.6] cm^g"\ which implies that the NFW-type 
cuspy profile persists in these haloes despite the SIDM in- 
teraction. 

Conclusions drawn from current N-body simulations of 
SIDM halo profiles may, therefore, require revision in light of 
our results. The scattering cross-section a ~ 218.5 cm'^g^^ 
is in an interesting regime which has not been studied be- 
fore by N-body simulations. This value may also help to re- 
solve the problem identified by the N-body simulations for 
smaller cross section, in which larger-mass haloes result in 
relatively larger cores. The small cross section regime corre- 
sponds to Q < Qth, for which the effect of SIDM conduction 
increases with increasing Q. As shown in Figure [TTI larger- 
mass haloes typically have larger Q-values, since their larger 
sizes more than offset the lower mean densities which result 
from their later formation. According to our results, how- 
ever, haloes from dwarf galaxies to clusters are not in the 
small cross section regime (i.e. low-Q regime). In fact, the 
dwarf galaxy rotation curves prefer Q = Qth, so large mass 
galaxies and clusters have Q > Qth, in general (see Fig. Illl 'l. 
According to Figure 

M this high-Q regime suppresses con- 
duction, yielding smaller cores (i.e. higher central densities) 
for higher-mass haloes. We predict, therefore, that as long as 
a ~ 218.5 cm^g~^ is used in an SIDM N-body simulation, a 
constant value, independent of velocity, will suffice to match 
both dwarf galaxy rotation curves and the mass profiles of 
larger-mass haloes. 



6 CONCLUSION/DISCUSSION 

CDM particles have been assumed to be coUisionless in the 
standard CDM theory of cosmic structure formation. De- 
spite the success of this standard theory, the elementary 
particle physics theory necessary to explain the origin and 



microscopic properties of the particles which comprise this 
dark matter is not yet known. It is natural for us to ask if the 
microscopic nature of CDM particles might lead to further 
constraints on this theory by astronomical observation. We 
have focused here on one such microscopic property, that of 
self-interaction by elastic scattering, and its effect on the in- 
ternal structure and dynamical evolution of virialized CDM 
haloes during galaxy and large-scale structure formation. 
The apparent discrepancy between the observed density pro- 
files of the haloes of dark-matter dominated dwarf and LSBs 
and those predicted by N-body simulations of coUisionless 
CDM may well be resolved if one assumes that CDM par- 
ticles interact with each other non-gravitationally. The self- 
interacting dark matter (SIDM) hypothesis is an attempt to 
produce such a nongravitational interaction between dark 
matter particles. 

We have derived the first fuUy-cosmological similarity 
solutions for halo formation in the presence of coUisionality. 
This provides an analytical theory of the effect of the self- 
interacting dark matter (SIDM) hypothesis on CDM halo 
density profiles as follows: 

• We have adopted the spherical infall model of cosmolog- 
ical halo formation, guided by the results of N-body simu- 
lations of CDM. The coUisional Boltzmann equation for the 
evolution of the gas of CDM particles yields a set of fluid-like 
conservation equations under the assumptions of spherical 
symmetry and isotropic velocity distribution. The effect of 
self-interaction collisions is accounted for by an effective con- 
ductivity term in the energy equation. This conductivity is 
valid for arbitrarily large or small collision mean free path 

Amfp. 

• For an Einstein-de Sitter universe (or a ffat universe 
with cosmological constant, at early times when matter 
dominates), the nonlinear growth of perturbations which 
leads to halo formation in the spherical infall model can be 
described by similarity solutions in the absence of conduc- 
tivity. In the presence of SIDM conductivity, self-similarity 
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is still possible, but only for mass perturbations 5M/M oc 
M~^^^ . Remarkably, this self-similarity required in our so- 
lution is well-motivated and justified by the theory of halo 
formation from peaks of the Gaussian random noise density 
fluctuations ( Hoffman fc Sh aham 1985 ). For galactic haloes, 
which form from density fluctuations whose power spectrum 
can be approximated by a power law P(k) oc the con- 

ditions required for this particular self-similarity naturally 
arise. 

• According to our similarity solutions, collisions of SIDM 
transport heat from the hotter, outer halo region into the 
colder core region. This process flattens the central density, 
and continuous infall pumps energy into the halo which sta- 
bilizes the core against gravothermal catastrophe. 

• These solutions are characterized by a single dimension- 
less quantity, the coUisionality parameter Q = api,rvir oc 
'"vir/Amfp, where a is the scattering cross section per unit 
mass, pb is the mean matter density, r^ir is halo virial radius 
and Amfp is the collision mean free path. The maximum flat- 
tening of central density occurs for an intermediate value of 
Q, Qth ~ 7.35 X 10"^, at which the halo is maximally re- 
laxed to isothermality. The density profile of the Qth solution 
matches that inferred from the observed rotation curves of 
dwarfs and low surface brightness galaxies (LSB) very well. 

• In the low-Q regime {Q < Qth), fiattening of the cen- 
tral density profile becomes stronger as Q increases. Pre- 
vious cosmological SIDM N-body s imulations with a ~ 
[0.1 - lOjcm^K" ^ lie in this regiin e iYoshida et al.ll2000bl 
iDave et alJboOlLIColin et alJl2002l) . Central density profiles 
became flatter as they increased the scattering cross-section, 
which is equivalent to increase in Q, because Q oa a. On the 
contrary, in the high-Q regime {Q > Qth), flattening of the 
central density becomes weaker as Q increases. This hap- 
pens because the scattering mean free path becomes shorter 
as Q increases. SIDM simulations which adopt a fully coUi- 
sional limit to derive the maximal de nsity flattening, which 
corresponds to ord inary gas dynamics jYoshida et alj2000al : 
iMoore et alj2000ll . report that haloes obtain density profiles 
with central cusps as steep as or steeper than those in col- 
lisionless N-body simulations. This seemingly puzzling be- 
haviour is easily explained: a —> oo corresponds to Q ^ oo, 
and therefore density flattening becomes negligible. 

• Under the assumption that dwarfs and LSBs formed 
at their typical collapse epoch in ACDM, a ~ 200cm^g~^ 
makes Q = Qth, much higher than previous estimates, a ~ 
[0.5 — 5]cm^g~^, based on N-body experiments. This value 
of cr, independent of halo mass, would make Q > Qth for 
clusters, which typically formed only recently, resulting in 
relatively less flattening of their central density profile and 
a smaller cor e. A velocity dependen t cross-section, ct oc l/n, 
suggested bv lYoshida et all i2000bri is thus unnecessary. 

• According to our similarity solutions, the solution for 
Q = Qth represents the solution most relaxed to isothermal- 
ity inside the virialized postshock region. It is notable, there- 
fore, that the Q = Qth solution is very similar to t he non- 
singular TIS solution of lShapiro. Iliev fc Ragal (^9^ for the 
post-collapse equilibrium structure of virialized haloes. The 
latter yields a mass profile almost indistinguishable from the 
mass profile of the Burkert fit to the rota tion curves of dwarf 
and LSB galaxies lllliev fc Shapirdr200 j) . as is the Q = Qth 
SIDM profile we have derived here. This suggests that the 
TIS halo model, which assumes that haloes are isothermal. 



is a natural outcome of the dynamical formation of CDM 
haloes when conductivity causes the halo to relax maximally 
toward isothermality. 

One may improve upon this work by performing a cos- 
mological N-body simulation. As our analysis is based upon 
self-similarity, or a constant mass accretion rate ^ -"-^ — g 

' a log a ' 

when the mass accretion rate deviates from this canoni- 
cal rate, our analysis is no longer valid. Several authors 
have investigated a realistic halo formation history both by 
an analytic approach and by N-body experiments, track- 
ing the history of the "most massive progenitor (MMP)". 
IWechsler et al.l (|20o3) performed a cosmological N-body 
simulatio n and tracked the growt h of MMP mass in a ACDM 
universe. iNusser fc ShethI lll99gt) calculated the growth of 
MMP for a power-law power spectrum and Ivan den BoschI 
(2002) calculated it for a CDM power spectrum, both us- 
ing the extended Press-Schechter theory. These studies show 
that mass accretion starts with a fast, rapid merger and ends 
with smooth, continuous accretion. This trend is clearly seen 
in equation (I66II . where the logarithmic accretion rate is de- 
creasing linearly with increasing scale factor a. The accretion 
rate obtained in this way is different from what is expected 
from HS. For n = —2.5 in the matter-dominated era, we find 
from equation ^ that = 4 or = 6. This fast 

^ » log t o log a 

accretion rate captures only the early mass accretion epoch 
given by equation (I66II . and therefore, one should not apply 
this rate to the later epoch when mass accretion becomes 
negligible. This is the moment where the evolution of SIDM 
haloes deviates from the self-similarity we assumed. 

For cluster-mass haloes, there is another reason that 
our assumption of self-similarity breaks down and should be 
improved in future work; the value of n which enters the 
self-similar infall model in ACDM actually depends weakly 
upon halo mass. The value we have adopted to ensure self- 
similarity in the presence of SIDM coUisionality, n — —2.5, 
is appropriate for the entire range of galactic halo masses. 
As Figure 121 shows, however, n increases with mass, and for 
clusters with mass above 10^" Mq, n > -2. If n / -2.5, 
self-similarity is violated by the presence of SIDM terms in 
the equations. Since the accretion rate is lower if n is higher 
(i.e. S^2sM_ ^ 3 N ^Yie flattening effect of SIDM on the 

^ a log a n + 3 ° 

halo central density proflle may be lower on cluster scales 
than expected from our self-similar solutions for n — —2.5. 

Does our prediction of high value of (t, a ~ 200cm^g~^, 
affect the abundance of dark matter substructure? Cravita- 
tional lensing flux anomalies have been interpreted as evi- 
dence for th e existence of da rk m atter substructures in the 
parent halo jMetcalf fc Z hao 2002: Dalai fc Kochanck 20021: 
iKeeton. Gaudi fc Pettersi„200a: .Mao et al 2004) . The self- 
interacting dark matter, if real, would suppress the number 
of dark matter substructures to some extent. We are not 
sure at this stage, however, how strong this effect will be: the 
simple assumption of cold, continuous mass infall prevents 
us from making any strong prediction. We instead describe, 
in the following discussion, the complexity relating to two 
main mechanisms for suppressing the substructure forma- 
tion when dark matter is coUisional: ram-pressure stripping 
and thermal evaporation. One should also note that the in- 
terpretation of the lensing flux anomalies is not settled yet. 
Recent analysis indicates that this anomaly may not require 
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Figure 13. Comparison of length scales for tlie Q = Q^^ case. 
The mean free path Amfp is comparable to the gravitational scale 
height H in the core region. However, the ratio of the mean free 
path to the gravitational scale height becomes larger as the radius 
increases. See also Figure |5] which includes other values of Q. 



any substructure in the primary len s halo, but onl y the sub- 
structure in the intergalactic space l)Metcaljl206^ V 

We now point out some caveats found in previous analy- 
ses, which also require more realistic, cosmological analysis. 
As many previous analytical estimates have been based on 
either isolated haloes or very simplified models, we believe 
that fully cosmological N-body simulations in the high-Q 
(a) regime should be carried out to clarify this issue. For 
instance, the restriction coming from t he susceptibility of 
SIDM haloes to ram-pressu re stripping (iFurlanetto fc Loebl 
120021: iNataraian et all2002tl may be relaxed if we remove the 
simplicity in their analysis. They determine the truncation 
radius of a galaxy in a cluster by the condition 



pc{r)vg = Pg{rt)c 



(79) 



where Vg is the velocity of the galaxy inside a cluster, Pc(r) is 
the density of the cluster at radius r, Pg{rt) is the density of 
the galaxy at its truncation radius rt, and ag is the internal 
velocity dispersion of the galaxy. Equation II79II is valid for 
highly coUisional fluid. They then use the restriction that 

Amfp ^ 



rt 



<7E(rt) 



< 1 



(80) 



where A^fp is the collision mean free path and S(rt), which 
is 0.024 g cm~'^ for their fiducial case, is the surface density of 
a galaxy at Vt- However, this is a crude way of describing the 
fluid regime where the equation (1791 can be applied. There 
exists a regime where the SIDM can be treated as coUisional 
only at the centre because the ratio of the mean free path to 
the gravitational scale height becomes larger as the radius 
of a galaxy increases (see Fig. I13| : when we express these 
length scales with dimensionless terms used for our similar- 
ity solutions, the scattering mean free path is Lmfp 

/3F/2, 



1 

DQ ' 



and the gravitational scale height is H 



In such 



cases, a "global" application of equation 17911 is not valid. 
Ram-pressure stripping in this case would not be as severe 
as in the case of purely coUisional fluid, because cluster 
SIDM particles have a high probability of penetrating an 
SIDM galactic halo deeply without experiencing collision at 
r > rt. Therefore, to apply equation I79II . we have to be more 



conservative in defining the fluid regime, which wi ll relax 
the constraint a<42cm'^g~^ iNataraian et al]l2002l) set by 
equation (EIB. 

iGnedin fc Ostrikerl (|2001j) constrained a from their nu- 
merical and semi-analytical calculation of the evaporation 
time of elliptical galaxies embedded in a cluster environ- 
ment. According to their analysis, such galactic haloes can 
evaporate within the age of the universe if tj = [0.3 — 
10'']cm^g~^. However, their analysis is based on a fixed clus- 
ter environment. Consider a cosmological variant of this 
problem: for instance, an elliptical galaxy may form by a 
recent merger at z = 1 when the temperature of the cluster, 
which only forms at 2 ~ 0, is very low. When the evap- 
oration time is about the age of the universe (~ 10^°yr), 
the elliptical galaxy has a fair chance to survive because the 
evaporation time is greater than the time from its forma- 
tion epoch to the present. The excluded range of a would 
then change substantially by shifting the marginal values, 
a ~ 0.3cm^g~^ or a ~ lO^cm^g"^, which correspond to 
the evaporation time of the order of the age of the universe. 
Moreover, their analysis is only valid for either the very long 
mean free path limit or the very short mean free path limit. 
They exclude the intermediate regime at ct ~ 200cm^g~^ 
simply by an extrapolation of these two regimes. 

Because of these problems, we assert that a more re- 
fined, fully cosmological analysis and new cosmological N- 
body simulations with a wider range of a values, including 
a ~ [100 — 500]cm^g~^, should be performed. Even though 
our analysis here is fully cosmological, it is restricted by the 
fact that it is based on a constant logarithmic mass accretion 
rate (i.e. = ^ = 7^ = GioTn = -2.5) that provides 

self-similarity. However, a more realistic mass accretion his- 
tory constructed from merger trees in N-body simulations 
shows a gradual decrease of the logarithmic mass accretion 
rate over time, as seen in equation 1661 (e.g. IWechsler et alJ 
2002). When the mass accretion rate becomes very small at 
late times, the underlying halo properties will deviate sig- 
nificantly from self-similarity. Moreover, the relatively high 
scattering cross-section which we find provides the best- 
fitting to dwarf galaxy rotation curves - ct ~ 200cm^g~^ 
- has never been tested in cosmological SIDM N-body sim- 
ulations. We will explore these issues further in future work. 
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APPENDIX 

We show how we got rioff for different mass scales. It is 
slightly different from the usual way to obtain UeS by differ- 
entiating the rms mass fluctuation, cta/. 
Typically, one has 



3 1-1- 



d In aM 
d\nM 



(81) 
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where 

Um-Mf) I r°° o , 

where m is a mass enclosed by a sphere of radius R which 
also defines the unperturbed mass M through 

M = ^R'po, (83) 

where po is the present matter density and the average { ) is 
taken over all positions of the centre of these spheres. This 
"top-hat" filtering results in a window function 



W{X) = — (sin(X) -Xcos(X)). 



(84) 



It is then straightforward to calculate riofF as a function of 
M using equation 18111 . 

However, we are interested in the rics which is valid if 
one considers the initial average overdensity around density 
peaks, Ao{R) (HS). Ao{R) is given by 



P{k)W{kR)k^dk, 



(85) 



where M and W{X) are defined by equation I83II and 18411 . 
respectively. From equation 1831 and equation 1851 . we can 
see that for a power- law power spectrum P{k) oc k", 

Ao{R) oc oc M-("+^)/^ (86) 

Therefore, one can obtain as follows: 



Weff 



-3 1 



dlnAo(A//) 
dlnM 



(87) 



The only difference in obtaining Ues is then the power 
of the window function W{X) in equations 1821 and 1851 . 
One can therefore expect that there would be only a minor 
difference, which we confirmed (see Fig. EP. We used the 
fitting formula for the transfer function bv lEisenstein fc HrJ 
il999tl . and worked in the ACDM concordance model with 

= 0.73, flmfi = 0.27, h = 0.7, erg = 0.9 and untilted 
Harrison-Zel'dovich primordial power spectrum. 
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